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Abstract 

We discuss the ill conditioning of the matrix for the discretised Poisson equa- 
tion in the small aspect ratio limit, and motivate this problem in the context 
of nonhydrostatic ocean modeUing. Efficient iterative solvers for the Poisson 
equation in small aspect ratio domains are crucial for the successful develop- 
ment of nonhydrostatic ocean models on unstructured meshes. We introduce 
a new multigrid prcconditioner for the Poisson problem which can be used 
with finite element discretisations on general unstructured meshes; this prc- 
conditioner is motivated by the fact that the Poisson problem has a condition 
number which is independent of aspect ratio when Dirichlet boundary con- 
ditions are imposed on the top surface of the domain. This leads to the first 
level in an algebraic multigrid solver (which can be extended by further con- 
ventional algebraic multigrid stages), and an additive smoother. We illustrate 
the method with numerical tests on unstructured meshes, which show that 
the prcconditioner makes a dramatic improvement on a more standard multi- 
grid prcconditioner approach, and also show that the additive smoother pro- 
duces better results than standard SOR smoothing. This new solver method 
makes it feasible to run nonhydrostatic unstructured mesh ocean models in 
small aspect ratio domains. 

Keywords: Finite elements, multigrid, unstructured meshes, small aspect 
ratio 
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1. Introduction 



There are many processes in tlie ocean (sucli as separating Western bound- 
ary currents, and density overflows) which are small in scale and restricted 
to particular regions, but which form crucial components of the global ocean 
circulation mechanism. It therefore seems attractive to design ocean models 
which use the flnite element method on fully unstructured meshes in order 
to incorporate some of these smaller scale features into a global ocean model 
(see Pain et al. (2005) for background and references). However, there are 
numerous pitfalls to negotiate in order to achieve this goal, arising from the 
fact that the global ocean is very thin: the horizontal lengthscale is thousands 
of times larger than the vertical lengthscale. 

One particular issue arises if one wishes to relax the hydrostatic approx- 
imation (Pedlosky, 1987), allowing for a model which is valid on both small 
and large scales. The nonhydrostatic pressure is obtained by solving a three 
dimensional elliptic problem with very large eigenvalues resulting from the 
horizontal scales, as well as very small eigenvalues resulting from the ver- 
tical scales. This means that the system is very ill conditioned. Since the 
Conjugate Gradient method, which is typically used for flnite element dis- 
cretisations with many degrees of freedom, has a convergence rate which 
scales with the square root of the condition number (see Shewchuk (1994) 
for example), this can have a catastrophic effect on the performance of the 
numerical model. 

In the ocean modelling context this problem was flrst encountered by 
Marshall et al. (1997). A solution strategy was proposed using a vertical 



preconditioner which solves the vertically integrated (aspect ratio indepen- 
dent) equations and then distributes the solution throughout the mesh. It 
was shown that the use of this preconditioner resulted in nonhydrostatic 
simulations which were as fast as hydrostatic simulations at the same resolu- 
tion. This strategy, has since been used in a number of nonhydrostatic ocean 
models, including those on horizontally unstructured grids such as [Fringer 



et al. (2006). However, the vertical averaging depends on the computational 



mesh being organised in vertical layers. This prohibits more general types of 
vertically unstructured meshes which may be required for multiscale simula- 
tions in which a small scale process is resolved within a large scale flow, or 
for hybrid meshes which accommodate both terrain following and isopycnal 
(constant density) layers. 

In this paper we extend the vertical averaging strategy so that it can 
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be applied to vertically unstructured meshes of large-scale ocean modelling 
such as those that can be used in the Imperial College Ocean Model (ICOM) 



(Piggott et al. , 2008). The extension is formulated by using the vertical ex- 
trapolation operator, which takes any point in the domain and returns the 
value of a function at the top surface directly above that point. This operator 
is the dual of the vertical integration operator, and can easily be approxi- 
mated on a vertically unstructured mesh. This extension is described within 
the context of the algebraic multigrid method; the framework also shows how 
to incorporate further algebraic multigrid stages into the "smoother" , which 
reconstructs the solution of the vertically integrated equations throughout 
the domain. This turns out to be necessary when a genuinely multiscale 
mesh is used in which the aspect ratio becomes 0{1) at the smallest scales. 
The aim is to obtain a numerical solver which has a convergence rate which 
is independent of the aspect ratio. 

The rest of this paper is organised as follows. In section [2| we describe 
the type of problems we wish to solve on small aspect ratio domains, and 
motivate them using the ocean modelling applications. In particular, this 
section explains why we cannot avoid solving an elliptic problem with Neu- 
mann boundary conditions on all surfaces, which is precisely the problem 
which gives rise to ill-conditioning. In section [3] we formulate the problem as 
a finite element approximation, in order to fix notation, and in section |4] we 
compute some estimates on the condition number for the Neumann bound- 
ary condition well as the case where Dirichlet boundary conditions 
are imposed on the top surface. It is observed that imposing the Dirich- 
let boundary conditions removes the small eigenvalues, and this motivates a 
preconditioning strategy in which one first eliminates the interior degrees of 
freedom to obtain an equation for the solution on the top surface, then one 
uses this surface solution as a Dirichlet boundary condition to reconstruct 
the solution throughout the domain. This paves the way for section [5] in 
which our proposed preconditioner is introduced, in the context of algebraic 
multigrid preconditioners for the Conjugate Gradient method. The precon- 
ditioner is tested in various examples in section [6j Finally, in section [7| we 
give a summary and outlook. 

2. Background: oceanographic applications 

In this section we describe how the pressure Poisson equation arises in 
nonhydrostatic models, which motivates the need to develop efficient solvers 
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for this equation in small aspect ratio domains. We shall also explain the 

types of boundary conditions that are imposed, in particular we shall explain 
why we need to tackle the problem of solving the pressure Poisson equation 
with Neumann boundary conditions. 

2.1. Nonhydrostatic equations 

The nonhydrostatic Euler-Boussinesq equations for a rotating stratified 
fluid on an /-plane are 

Po +U-VU + fz xi?j = -Vp - gpz, (1) 

V-u = 0, (2) 

dT 

— + u-vr = 0, (3) 

where u is the velocity, po is the (constant) reference density, z is the unit 
vector in the upward direction, / is the Coriolis parameter, p is the pres- 
sure, g is the gravitation constant, T is the potential temperature and p is 
a prescribed function of the temperature. Here we do not complicate the 
exposition by including viscosity or diffusivity terms, or dependence of the 
density on salinity or pressure. 

We solve the problem in an "ocean shaped" domain Q with bottom bound- 
ary dflfiooT, coastal boundaries Sficoast and top boundary dfltop (which may 
be allowed to move up and down to accommodate surface waves). The hori- 
zontal extent of the domain is L, and the vertical extent of the domain is H; 
in this paper we concentrate on the difficulties when the aspect ratio H/L 
of the domain is very small (H/L Ki 1/1000 for an ocean basin). We shall 
parameterise the top surface by 

z = r)(x,y), 

with 7] = when the fluid is at rest, and make the additional simplifying 
assumption that the coastal boundaries are vertical. 
We consider two types of boundary conditions: 

• Rigid lid: 

u-n = 0, X e dilfiooT U SOcoast U dUtop, 
V = 0, 

for a constant top surface height z — 0. 
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• Free surface: 

U-n = 0, X e Sfifloor U dVtcos.su 
P = Pa, X e dQtop, 

r]t = -Uh ■ \/hV -w = -. (4) 

n ■ z 

In the subsequent section we shall show that both of these equations effec- 
tively result in an ill conditioned pressure equation in the small aspect ratio 
limit; this always occurs for the rigid lid case and also occurs for the free sur- 
face case when one wishes to take large timesteps and has an unstructured 
mesh. 



2.2. Solving for the pressure 

The role of the pressure in these equations is as a Lagrange multiplier 
which enforces the incompressibility condition (|2| (in fact it enforces that 
the solution to + V ■ (Du) = is D = 1, and equation ^ is then a direct 
Holm et al. ( |1998 )) for example. To solve for the pressure, 



see 



consequence 

take the divergence of equation (jlj): 



V ■ (po (it ■ V + fzx)u + zgp) 



which is a Poisson equation for the pressure p, given the other variables u and 
T. In practise, one usually solves for the pressure update po(l>/^t = p""*"^ — 
(where is the pressure at time level n) using one of the families of pro- 
jection methods based on Chorin ( 1967[ ); Temam (1969). These methods 
are typically predictor-corrector schemes in which a predictor u* is obtained 
using the pressure from the previous timestep, without enforcing the incom- 
pressibility condition ([2]), and then constructing a correction 



u 



n+l 



U 



V0, 



(5) 



subject to 

V ■ = 0. 
Taking the divergence of equation ^ gives 



-V u*. 



(6) 



For more details, see Karniadakis and Sherwin (2005); Gresho and Sani 
(|2000l). 
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On slip boundaries where u n = 0, since u* already satisfies the boundary 
conditions, then 

so that it""^^ keeps the same boundary conditions. In the free surface case, 
the pressure must stay constant on the surface, so satisfies the zero Dirichlet 
condition on the top surface: 

— = on Silcoast U dQaoor, = OU dQtop- 

However, the equations support free surface (barotropic) waves which are 
very fast compared to the other dynamics. Moreover, the highest frequency 
waves have small amplitudes and we are not concerned with the details of 
their evolution in large scale models. The standard ocean modelling approach 
is to apply a splitting method in which one integrates the barotropic waves 
in time using a much smaller timestep; the remaining (baroclinic) dynamics 



is then integrated using a larger timestep (see Shchepetkin and McWilliams 



(2005 ) for example). However, if the mesh used is unstructured in the vertical 
direction [i.e. the mesh is not arranged into layers) then it is not possible to 
obtain such a splitting. An alternative approach which is often used in coastal 



engineering applications (see Labeur and Pietrzak (2005) for example), is to 



construct a new "piezometric" pressure 

Pix, y, z) = p{x, y, z) - pQg^{x, y) - Pa, 
which specifies the Dirichlet boundary condition 

P = -Pogv, on dQtop- (8) 

This means that we do not need a separate free surface variable as the value 
can simply be read from p at the free surface. This piezometric variable 
satisfies the same pressure Poisson equation with modified right-hand side 
and different boundary conditions. 

Taking the time derivative of equation ([s]) and substituting the kinematic 
boundary condition ^ gives 

dp un 
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If we wish to take large timesteps compared to the timescale of the free surface 
waves, the barotropic-barochnic spht is not available on unstructured meshes 
and it is necessary to solve the equations using a linearly implicit method. 
Here we will describe the simplest case, the backward Euler scheme, but the 
development of higher order schemes is similar. The backward Euler scheme 
gives 

Substituting the pressure update gives 



g-^, on dntop, 



At2 ^ n ■ fc 

which is a Robin boundary condition for (p. The ratio of these two terms is 
approximately 



At2 



H f H 



9^ 



gAt^ \cAt 



where c is the barotropic wave speed y/gH. This is the square of the ratio of 
the time it takes a barotropic wave to travel a distance H to the timestep; 
if we wish to take large timesteps then this quantity is small and we recover 
the Neumann boundary condition ([T]). 

All of this means that if we wish to take large timesteps with an unstruc- 
tured mesh, then we must have an efficient method for solving equation ^ 
with boundary conditions ([T]). In this paper we shall see that this equation is 
ill conditioned when the aspect ratio H/ L is small. We will introduce a new 
multigrid preconditioner which allows efficient iterative methods for solving 
numerical discretisations on this problem on unstructured meshes. 

2.3. Hydrostatic equations 

In contrast, if one makes the hydrostatic approximation which is used in 
many ocean models, 

Pz = -gp, (9) 

then the pressure can be obtained by integrating this equation, using the 
boundary condition p = on the top surface dQtop- On a vertically struc- 
tured mesh this equation can be accurately integrated from top to bottom in 
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columns. On an unstructured mesh, one approach is to differentiate equation 
^ to obtain the eUiptic problem 

O^p p Op 

-r^ = -9q^, P = on d^top, -r^ = -gp on dnnoor- (10) 



The hydrostatic pressure equation ( 10 ) does not suffer from the same ill con 



ditioning as the nonhydrostatic pressure equation ^ with boundary con- 
ditions ([T]), since the smallest eigenvalue is independent of the aspect ratio 
e (as we shall see later). This provides the benchmark for nonhydrostatic 
pressure equation solver methods, with the aim of developing methods for 
the nonhydrostatic equation which are as fast as methods for the hydrostatic 
equation. 



3. Finite element formulation 

In this paper we consider the finite element approximation to the Poisson 
equation, 

-V^<^ = /, 

on the domain Q with homogeneous Neumann boundary conditions 

on 

We consider a domain Q with horizontal scale L, vertical scale H, horizontal 
velocity scale U, and choose nondimensional domain coordinates 

x' = x/L, y' = y/L, z' = z/H, 

In these coordinates, the Poisson equation becomes 

with boundary conditions 

(e^n;, e^n'a, n'g) ■ V'0 = 0, on dVL\ n' = {n[, n'^, n'3) (12) 

where e = H/L is the aspect ratio of the unrescaled domain Q, and fl' is the 
rescaled domain with normal n'. Henceforth we drop all the primes. 
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Given a solution of equation (11) subject to boundary conditions (12), 
it is possible to obtain a whole family of solutions 



+ c, c G 



and hence the solution is not unique. However, the constant c does not have a 
physical effect since the pressure only appears in the Navier Stokes equations 
as a gradient. This means that we may arbitrarily fix this constant. This is 
usually done by requiring that 



or 



.d\/ = 0, 



M = 0, 



(13) 



(14) 



for some chosen point a^o G fi. Condition (14) is easier to implement in the 



finite element method on a general unstructured mesh, and condition (13) 
can subsequently be imposed by subtracting off the integral of (p. Hence, we 
shall require condition (14). Here we shall additionally require that Xq G 



Sfitop! the approach we describe does not depend on this requirement but it 
simplifies the exposition. It was noted in Bochev and Lehoucq (2005) that 



good performance can also be obtained with the Conjugate Gradient method 
if the value of the constant mode is not fixed in the assembled equations, but 
instead projected out each iteration of the CG solver, but we do not discuss 
that case in this paper. 

A finite element discretisation of this equation is obtained by writing the 
weak form of equations (11) with boundary conditions (12). First we define 
the function space 

H\xo) = {(f)eH': (j){xo) = 0}, 
where is the Sobolev space with norm 

Hrm= [ \<p\'+m'dv, 

Jn 

and we seek G H^{xq) such that 

5,(7^,0) = F(^), (15) 
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for all test functions G H^{xq), where 



n dz dz 



dip d(f) dip dcj) 
dx dx dy dy 



dV 



To construct the Galerkin finite element discretisation of these equations 



(see Brenner and Scott (1994) for example), we select a finite dimensional 
trial space V{xq) C H^{xq) (typically by constructing a polygonal mesh on 
the domain Q and constructing piecewise polynomials), and seek (p^ G ^(a^o) 
such that 

B,iiP',(P') = Fiip'), (16) 

for all test functions ip^ G V. If we instead wish to solve the equations with 
a Dirichlet boundary condition imposed on the top surface dfltop, 



g\ yx G on 



top; 



(17) 



then we construct the homogeneous Dirichlet space V out of functions in V 
which are zero, 



v{dn 



top) 



G V, (i)\x) = yxe oritop} 



and seek (p G V{dQtop) such that 



B.{ip\x') 



(18) 



for all test functions ip G V(9fitop), where is some chosen function which 
satisfies (17). The total solution is then cp^ = cp^ + ■ 



This is often 



.mine referred to as "lifting" the boundary conditions; see Karniadakis and 



Sherwin (2005, for example) for more details. We note that this becomes the 



referred to as "lifting" the boundary conditions; see Karniadakis 



and Sherwin (2005) for example for more details. We note that this becomes 



the lllllll .r779 Galerkin finite element approximation of the hydrostatic 
problem (10) with suitably chosen and when e = 0. 

To solve these equations we expand the trial functions 0^, the test func- 
tions ip^ , and the right hand side function (or — x^ in the surface 
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Dirichlet case) in basis function expansions for V (or V in the surface Dirich- 
let case) 



X 



i=l 



i=l 



i=l 



and substitute into equation ( 16 ) (or equation ( 18 )) to obtain a matrix vector 
equation 

A,cf> = Mf, (19) 

where 



A 



dN, dN, 



dz dz 



dNi dNj dNi dNj 
dx dx dy dy 



d V, Mr 



NiNjdV. 



Here is a positive definite sparse matrix, and we wish to solve equation 
(19) iteratively using the preconditioned Conjugate Gradient method. In 
this paper we shall choose a basis expansion of V{xq) so that the vector (f) 
of basis coefficients of d) takes the form 







(20) 



where (p G R"* is the vector of coefficients corresponding to basis functions 
which are zero on dQtop, and 0' G is the vector of the coefficients cor- 
responding to the remaining basis functions. These are typically chosen to 
vanish on every finite element "node" which is not on dQtop. In this ordering, 
we write 



CT A, 



Mf 



(21) 



Note that the matrix obtained from equation (18) (which we denote as A^), 
is a sub-matrix of the matrix obtained from equation ( 16 ) (which we denote 
as A^). This matrix solves the problem with Dirichlet boundary condition at 
the top surface, rather than Neumann. We shall make repeated use of this 
decomposition throughout the rest of the paper. 



4. Condition number estimates 

In this section we obtain estimates in the small aspect-ratio limit for the 
condition number of the matrix A^ which we developed in the previous sec- 
tion, for both the Neumann and Dirichlet boundary condition cases. When 
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the condition number is large, the iterative method can be very slow to con- 
verge, and so it is important to understand the dependence of the condition 
number on the aspect-ratio. We shall note that the Neumann boundary 
condition case (which is the case of interest for oceanographic problems) has 
a condition number which scales like e~^, whereas the Dirichlet boundary 
condition case has a condition number which is independent of e as e — > 0. 
This motivates our proposed preconditioner. 

In this section we estimate the minimum and maximum eigenvalues of 
symmetric matrices by using the Rayleigh quotient estimates 

>^min — mill 7^ , ArvioY — IHRX 7^ . 

To facilitate these estimates, we define the vertical operator Aq and horizontal 
operator Ah with coefficients 

_fdN,dN 

^H,ij 



dx dx dy dy 



so that 

A, = Ao + e^An. 

We first construct an upper bound for the minimum eigenvalue \inhi of 
A^. First we note that the intersection of the null space of Aq with Ah is the 
zero vector. If y G ker(74o) and y e 'kex{AH), then y e 'k&!:{A^). However, 
is invertible so t/ = 0. This allows us to estimate the minimum eigenvalue: 

_ . <t>^A,ct> 

(tP'A,<t> 



<i min 

^ 2 • 4>^AH<t> 

< e mm =^ , 

<^^0,<^eker(Ao) (f) (f) 



where 



ct>^AHcf> 

Co = mm , 

<^^0,<^eker(Ao) (f)^ (f) 
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which is bounded away from zero since ker(Ao) fl ker(A//) = {0}. Next we 
estimate the maximum eigenvalue Amax of A^. 

= max , 

> max — = ci, 

Here, ci is the maximum eigenvalue of Aq which is independent of e. Next 
we compute the condition number of A^ which is the ratio of the largest and 
smallest eigenvalues 

Cond(A,) = ^ > 

Amin Co 

This means that the condition number is unbounded in the small aspect ratio 
limit e ^ 0. Since the convergence rate of the Conjugate Gradient method 
typically scales with the square root of the condition number, this means 
that the Conjugate Gradient method becomes very slow as e — 0, and we 
must find a preconditioner which makes the condition number independent 
of e. 

In figures [T] and [2] we illustrate these estimates. Two tetrahedral meshes 
in a box domain were generated, one which is arranged in horizontal layers, 
and one which is fully unstructured in all three dimensions. The finite el- 
ement approximation to the Laplace equation was applied to these meshes, 
having rescaled the coordinates to various different aspect ratios. The eigen- 
values were then numerically computed using Arnoldi iteration. Figure [T] 
shows the eigenvalues for the layered mesh, with various different aspect 
ratios. We observe a gap in the spectrum between the eigenvalues corre- 
sponding to 2;-independent eigenvectors (we call these horizontal modes) and 
the eigenvalues corresponding to 2;-dependent eigenvectors. The size of this 
gap is proportional to e~^. It can also be observed that the ratio between 
the largest and smallest eigenvalues is proportional to e~^. In the unstruc- 
tured mesh, the distinction between the horizontal modes and the rest of 
the eigenvectors is less clear. The lack of horizontal alignment in the mesh 
means that there are numerical errors in the finite element approximation 
of the vertical derivatives which scale with the horizontal widths Ax of the 
elements fin this case Ax^ since we have used linear finite elements) and 



13 




Figure 1: Plot showing eigenvalues for the matrix A arising from the discretisation of the 
Poisson equation in a box domain with Neumann boundary conditions on all sides. The 
box was decomposed into a tetrahedral mesh divided into a number of horizontal layers 
(i.e. the mesh is structured in the vertical) and rescaled into various different aspect ratios 
e. Note that there are a cluster of small eigenvalues which are independent of e: these 
eigenvalues correspond to the z-independent eigenmodes. As e decreases to zero, the width 
of the spectral gap between these and the remaining eigenvalues scales in proportion to 



hence are the same order of magnitude (or larger) than the exact eigenvalues 
of the horizontal modes. Figure [T] shows the eigenvalues for the unstructured 
mesh with various different aspect ratios. We observe the same scaling 
for the ratio between the largest and smallest eigenvalues, but there is no 
spectral gap since the small eigenvalues are polluted by numerical errors in 
the vertical derivatives. 

In contrast, we note very different scaling behaviour when the Neumann 
boundary condition on the upper surface (rigid lid) is replaced by a Dirichlet 
boundary condition, resulting in the matrix A^. For this case, we shall obtain 
an upper bound on the condition number. First we bound the minimum 
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Figure 2: Plot showing eigenvalues for the matrix A arising from the discretisation of 
the Poisson equation in a box domain with Neumann boundary conditions on all sides. 
The box was decomposed into a tetrahedral mesh which is completely unstructured in the 
vertical and rescaled into various different aspect ratios e. Note that there is no longer a 
spectral gap, which is replaced by a more evenly spaced spectrum, but that the spread of 
eigenvalues still increases in proportion to e~^. 
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eigenvalue Amm of from below. 

— mm ^ — ^Jf^ — ) 

> mm , 

= C2, 

which is the minimum eigenvahic of and is bounded away from zero since 
^0 is non-singular. Next we bound the maximum eigenvalue Amax of A^. 

- _ 

'^max — max rp , 

<t>¥^o (f) 

0^ (Aq + e^An) ^ 
< 2 inax , 

<t>^0 (f) (f) 
= 2C3, 

provided that e is sufficiently small; here C3 is the maximum eigenvalue of 
Aq. This means that the condition number of ^4^ is bounded by 

Cond(/l,) = h^<'^ = 2Cond(Ao), 

Amin C2 

which is twice the condition number of Aq, and is independent of e (this is 
not a sharp estimate, but illustrates the scaling with e). 

The contrast between the condition number scaling of A^ and A^ moti- 
vates a preconditioner strategy in which one solves a reduced problem for 
the solution on the surface dfltop, and then uses this surface solution as a 
Dirichlet boundary condition to reconstruct a solution throughout Q. This 
reconstruction step amounts to inverting A^ which, as we have just seen, 
has a condition number which is independent of e. We shall describe this 
preconditioner strategy in the following section. 
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5. New preconditioner 



In this section we develop our new preconditioner for equation ( 19 ), which 
is derived from the strategy of ehminating the degrees of freedom associated 
with the solution in the interior of Q to give a reduced equation on the sur- 
face dQtop- In each iteration the preconditioner will approximately solve this 
problem, and then approximately reconstruct the solution in the interior us- 



ing A^. We describe the preconditioner as follows: in Section |5.1| we briefly 
summarise the general algebraic multigrid preconditioning strategy. In Sec- 
tion 5.2 we introduce a reformulation of equation the (19) that decomposes 



the inverse of A into a vertically lumped system and a system with a Dirichlet 



boundary condition on top. Sections 5.3 and 5.4 explain the approximations 



that need to be made to this decomposition to apply it as a preconditioner. 

5.1. Algebraic multigrid preconditioners 

The general idea of multigrid methods is to tackle multiscale, ill con- 
ditioned problems by trying to solve for the different components of the 
solution, associated with different length scales, separately. This is accom- 
plished by a sequence of coarsening operations, in which the dimension of 
the problem is reduced step by step. The coarser system no longer supports 
the smaller scale features and has therefore an improved condition number. 
Thus the large scale, small eigenvalue modes can be efficiently solved on a 
reduced system, whereas the small scale, large eigenvalue modes are easily 
reduced with standard preconditioners such as SOR (therefore in this context 
referred to as smoothers) . 

Classical geometric multigrid methods, implement this coarsening step 
via a coarsening of the mesh on which the problem is defined, for instance 
via a /i — >■ 2/i coarsening on structured meshes. Algebraic multigrid (AMG) 



methods (see Stiiben (2001 ) for an introduction), use the algebraic properties, 
matrix graph and coefficients, of the matrix to construct a coarsening oper- 
ator. This more general approach has the advantage that it works equally 
well for unstructured mesh discretisations. Additionally it is possible to take 
anisotropies in the problem into account by selecting only matrix graph con- 
nections associated with large matrix coefficients. 

For symmetric problems, to keep the problem symmetric at each level, 
the prolongation operator P, that maps the solution of the reduced system 
back to the previous level, is usually the transpose of the coarsening operator 
P^. Given the original matrix A at the finest level, the matrix of the reduced 



17 



system, is given by 

P'^APy = P'^b, y e M™, (22) 

After solving at the coarsest level, the solution is mapped back using y — >■ 
X = Py. 

The smaller scale modes that are only represented in the full problem are 
reduced by applying a smoother S. As this smoothing step will also again 
change the solution at the coarse level, i.e. the solve at the coarse level and 
the smoothing for the full problem are not independent, the whole procedure 
of restriction, coarse solve, prolongation, and smoothing needs to be applied 
in an iterative manner. To keep things symmetric, the smoothing step 5* 
after prolongation is usually mirrored by a transpose smoother S"^ before the 
reduction. For instance a forward sweep of SOR before the restriction can 
be accompanied by a backward sweep after the prolongation. 

A typical 2-level multigrid cycle (V-cycle) then looks like 




Finally by replacing the coarse solve with a multigrid V-cycle applied to 
the reduced system, the multigrid method can be extended recursively to 
multiple levels. 

Best results are obtained if the multigrid V-cycle is embedded, as a pre- 
conditioner, in a Krylov subspace method. Different multigrid precondition- 
ing approaches are formed by different coarsening strategies and different 
choices of smoothers. The algebraic multigrid preconditioner used in the re- 



sults section, implements the smoothed aggregation approach of Vanek et al 



(1996). This method is known to work very well for strongly anisotropic 
elliptic problems. As will be shown in the results section however the con- 
vergence rate is not independent of the aspect ratio. Therefore, we seek to 
improve upon this purely algebraic method in the following sections. 

5.2. Schur complement equation 

Motivated by the analysis in section |4| where we observed that the con- 
dition number of the linear system becomes independent of the aspect ratio 
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if the Neumann boundary condition on top is replaced by a Dirichlet condi- 
tion, we proceed by constructing a reduced system where we first solve for 
the solution 0' on dfl^op, and then reconstruct (cf. the decomposition of 
(f) into (f)' and 4> in (20) and (21)). This can be done by solving the Schur 

(23) 



complement equation 



{B - CA^ 'C^) 4>' = h' 

V ' 

Schur matrix 



and then solving 

A4, = -C^cj>' + b', (24) 

to reconstruct the solution in the interior. Note that dim(0') ^ dim(0), and 
also that the reconstruction equation (24) has a condition number which is 
independent of e as e — 0. 

The Schur complement matrix contains and hence is a full matrix 
which is expensive to assemble and solve. Hence we shall propose a strategy 
to form approximations to equations (23) and (24) which can be used as a 
preconditioner in a manner similar to a multigrid preconditioner. 

5.3. Extrapolation operator 

To build the approximation to the Schur matrix, we first note that equa- 
tion 23 may be rewritten as 



E^A,E(f)' = E^b, 



(25) 



where 



E 



We call E the extrapolation operator. Given 0' G 

= Ecf>', 



the operation 



produces the finite element discretisation of the solution of the Laplace equa- 
tion 

= 0, (26) 

with Neumann boundary conditions dip/dn = on the coasts and bottom 
surface, and Dirichlet boundary conditions 



(27) 
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on the top surface, where is the function on f2 with finite element basis 
function coefficient vector (p and 0' is the function on dfltop with finite ele- 
ment basis function coefficient vector 0'. Note that in the small aspect ratio 



limit, (26) converges to 



— , 



0, 



with (j) = (f)' on the top surface, and d(f)/dz = on the bottom surface. The 
solution is then the vertical extrapolation of 0' i. e. 

0(a;,2/, z) = (t)\x,y). 

We shall use this in subsequent sections to construct an approximation to E. 



We can eliminate using equation (|24|) to obtain 

= 



0' 



4 ' (-c^0' + b') 



Ecf)' 







that is 



A: 



E{E^A,E)-^E^b + 
E{E^A,E)-^E^ + 



A 



-1 



a: 



(28) 



smoother 



7 



The two parts of this formula may be interpreted in the context of a general 
multigrid strategy. The first term is similar to a 2-level multigrid cycle with 
prolongation operator E. It projects the equation to a vertically lumped 
system. The second term acts on the vertical modes in the solution and can 



therefore be seen as an additive smoother. It is to be noted that (28) provides 



an exact solution for the inverse of A^, provided both inversions are performed 
exactly. However, both E and A^ are dense matrices; in the next two 
sections we will provide approximations for both the extrapolation operator 



equation (19). 



E and the inverse of A^, such that (28) can be used as a preconditioner for 
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Figure 3: Vertical extrapolation performed by projecting each node in the full mesh (here 
depicted in 2D) straight upward onto the top surface mesh (here ID), and interpolating 
the value in this projected node from the surrounding nodes in the surface element. 



5.4- Approximation of the vertical extrapolation operator E 

As noted in the previous section, the operator E given by discretisation 
of (26) with Dirichlet boundary conditions condition (27) on top, converges 
to a vertical extrapolation operator as e — > 0. We therefore expect that for 
large e, such a vertical extrapolation operator, between the top surface mesh 
and the full mesh, is a good approximation of E. This operator can simply be 
constructed by projecting nodes of the full mesh in the vertical direction onto 
the surface mesh, and interpolating within the surface triangle each projected 
node lies (see figure [S). This gives an approximation E : R*"' — > R*" of E with 
a limited stencil: for a continuous linear (PI) discretisation in 3 dimensions 
it connects every interior node with three nodes of the surface mesh. 

Another approach to find a sparse approximation of E is given by project- 
ing E onto a chosen sparsity pattern using a mod i ficatio n of the symmetric 



sparse approximate inverse (SSPAI) (Benzi et al. 
writing 

E- 



1996). E is obtained by 
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and then selecting a sparsity pattern for F. The non sparse entries of F are 
then obtained by minimising 



\\F^A-cr, 

subject to the sparsity constraints, where || ■ || is the Frobenius norm. This 
leads to decoupled sparse matrix problems 

AiVi = Vi, i = 1, . . . , 71, 

where Vi is the i-th column of F restricted to the sparsity pattern of that 
column, Vi is the i-th column of C'^ restricted to the sparsity pattern, and 
Ai is restricted to the sparsity pattern of the i-th column. We do not 
pursue this approach in this paper, preferring to use the projection approach 
described above. 

5.5. Additive smoother 

By replacing E with E, we have produced the approximate inverse 

a;' ^ E{E^AE)-'E^ + a:' ^ . (29) 

To use this as a preconditioner, we must also approximate the additive 
smoother. The second term could be evaluated exactly by solving a ma- 
trix equation A^cp = b for the interior part of the residual. Although, as 
noted before, this system is much better conditioned than the full system, 
the solution of an elliptic equation on the interior of the mesh is still quite an 
expensive operation that needs to be performed during each application of 
the preconditioner within each the Krylov iteration. Moreover the solution 
of this interior equation needs to be done using an iterative Krylov method 
as well. It is well known that embedding a Krylov method within another 
Krylov iteration, requires the use of a flexible Krylov method for the outer it- 



eration {e.g. FGMRES (Saad, 1993)). A major drawback would therefore be 
that this approach would inhibit the use of the Conjugate Gradient method 
for the outer iteration. 

A very simple smoothing strategy is obtained by realising that the first 
term of the proposed preconditioner is just the first stage of a general multi- 
grid method. In this projection the long scale, horizontal modes are separated 
out and the vertical projection can therefore be interpreted as a general coars- 
ening step such as those in any multigrid method. The necessary smoothing 
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step to filter out the short scale modes, is there often done by application 
of one or more SOR iterations of the entire system. For small aspect ratio 
problems this may therefore be enough to reduce the vertical modes in the 
error. 

In some cases, the mesh may contain a lot of structure in the vertical as 
well as in the horizontal. For instance an adaptive mesh model might focus 
resolution on physics related to the baroclinic modes of the system. In such 
cases the simple SOR smoothing may not be enough. The vertical lumping 
step would take out too much structure in one step. This may be compared 
to so called "aggressive coarsening" techniques in general multigrid methods 
that are usually accompanied with improved smoothing techniques. A more 
accurate approximation of the second term in (29) would be to replace the 
inverse matrix by a full cycle of the general AMG method applied to 
A^. The next section will provide a comparison of the simple SOR smoother 
with this more advanced additive smoother. 



6. Numerical experiments 

In this section we present numerical results which test out our precondi- 
tioner on matrices obtained from the linear finite element approximation of 
the Laplace equation with the horizontal coordinates rescaled to various as- 
pect ratios. The solvers were developed using the open source PETSc library 



(Balay et al. 1997) 



To compute errors, we selected a right hand side for the matrix vector 
equation by choosing a solution and multiplying it by the matrix. This allows 
us to compute errors exactly at each iteration of the solvers. Throughout 
this section we use the inf-norm to measure the magnitude of the error: our 
rationale for this is that we are motivated by multiscale applications in which 
one may be very concerned with the numerical solution in one small region 
of the domain (for example one may wish to embed a convection cell in an 
ocean basin and observe how it is affected by the large scale dynamics). In 
this case it may be possible to obtain a small L2 error whilst the solution in 
the small region is still inaccurate. We present plots of error against number 
of iterations, and also against floating point operations (flops). The flop 
count is provided by a intrinsic PETSc routine. 

In this section we obtain results from two meshes, both of a 1x1x1 cube; 
the coordinates of the meshes are then rescaled to a range of aspect ratios 
in the small aspect ratio limit. Mesh A is a Delaunay triangulation for a set 
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of 57453 roughly equispaced points in the cube; the points are not arranged 
in layers and hence the mesh is unstructured in the vertical. Mesh B is 
a Delaunay triangulation on a mesh in which the majority of the points 
are clustered at the centre of the cube (there are 99017 points in total), 
leading to very small elements. This is a truly multiscale mesh in which 
small eigenvalues exist due to both the small aspect ratio and also due to 
small elements. Our aim is to develop robust, efficient matrix solvers for 
these challenging multiscale meshes. 

The three preconditioners that are compared are: 

• A general AMG method based on the smoothed aggregation method 



(Vanek et al. , 1996). This uses our own implementation constructed 
using the "MG" interface provided by PETSc. The smoothing at each 
level is given by a single forward SOR sweep {uj = 1.0) as a pre- 
smoother and a backward sweep for post-smoothing. The coarsening 
strategy is based on the strongly-coupled connection criterion 



\Aij\ > e AiiAj^ 



n 



where a e of 0.01 has been chosen. The smoothing in the aggregation 
operator uses oj = 2/3. 

The preconditioner given by the vertically lumped approach 



AJ^ ^ E{E^A,E)'^E^ 

where the vertically lumped system is approximately solved using a 
single multigrid cycle applied to E'^A^E. This is combined with a 
single forward and backward SOR sweep as respectively a pre and post 
smoothing step. Thus the vertical lumping operator E is treated as an 
ordinary coarsening operator, and the vertical lumping of the equation 
is simply the first of a multilevel multigrid cycle. 

• As a last approach, the above multigrid cycle, including the vertical 
lumping as the first coarsening step, is combined with the additive 
smoother, where A^ is approximated applying a cycle of the smoothed 
aggregation AMG method to A^. 

In all cases the full multigrid cycle is applied as a preconditioner within each 
iteration of the Conjugate Gradient method. 
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6.1. General multigrid methods 

In figure |4| we compare the smoothed aggregation preconditioner to two 
other algebraic multigrid algorithms available through the PETSc library 
(BoomerAMG and Prometheus), and the classical Symmetric Successive 
Over-Relaxation (SSOR) preconditioner. The four preconditioners are com- 
bined with the Conjugate Gradient method in solving the pressure Poisson 
equation on Mesh A rescaled to an aspect ratio of 1/1000 (a reasonable as- 
pect ratio for large scale oceanographic applications). It is shown that the 
smoothed aggregation is substantially more effective than the other meth- 
ods (we additionally notice that the other multigrid methods are not much 
more effective than SSOR in this small aspect ratio example, although both 
methods contain a number of parameters and it may be possible to obtain 
better results by tuning). All four methods produce long "plateaus" in the 
error that are maintained for many iterations before the error finally drops. 
This means that none of these preconditioners result in feasible methods for 
solving the pressure Poisson equation in small aspect ratio domains. 

In figure |5} we further illustrate the problems that occur when the Con- 
jugate Gradient method with the smoothed aggregation multigrid precon- 
ditioner is applied to the Poisson equation on mesh A. The error is plot- 
ted against the iteration number for the preconditioned Conjugate Gradi- 
ent method for various aspect ratios. As the aspect ratio of the domain 
decreases, the condition number increases and the convergence rate of the 
iterative solver gets slower. We observe a "plateau" in the error which is 
maintained for many iterations before the error finally drops; this plateau 
becomes longer and longer as the aspect ratio decreases. For small aspect 
ratios this makes the iterative solver prohibitively slow. 

6.2. Preconditioner with vertical lumping 

In figure |6| the same information (error plotted against iteration number 
for the solution of the matrix system obtained from Mesh A) is given for 
the vertically lumped preconditioner using an SOR smoother. We note that 
in contrast to the standard multigrid preconditioners tested in the previous 
subsection, there is no plateau and the convergence rate becomes indepen- 
dent of e for small aspect ratios. We attribute this fast convergence to the 
removal of small eigenvalues in (nearly) vertically-independent eigenmodes 
by the vertically lumped preconditioner. For small aspect ratios there is an 
exponential decay of error with iteration from the very first iterations. 
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Figure 4: Plot showing error versus computation time for the smoothed aggregation 
multigrid preconditioner combined with the Conjugate Gradient method, compared with 
two other multigrid preconditioners available within PETSc (namely BoomerAMG and 
Prometheus), and the SSOR preconditioner. The smoothed aggregation method uses a 
coarsening strategy which is weighted by the matrix entries. When the aspect ratio is 
small (it is 1/1000 in this case), this method tends to aggregate degrees of freedom which 
are nearly vertically aligned. The smoothed aggregation method is substantially more 
effective than the other multigrid methods and the SSOR method, but we still observe a 
long period during which the error is not reduced. 
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Figure 5: Plot showing error (using the c»-norm) against number of iterations, for the 
Conjugate Gradient method applied to the Poisson equation discretised on mesh A, using 
the smoothed aggregation multigrid preconditioner. The mesh has been rescaled to various 
different aspect ratios e as indicated in the plot. The number of iterations required to 
converge increases with decreasing e, with a long "plateau" for small aspect ratios. 
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Figure 6: Plot showing error (using the cxi-norm) against number of iterations, for the 
Conjugate Gradient method apphed to the Poisson equation discretised on mesh A, using 
the vertically lumped preconditioner. The mesh has been rescaled to various different 
aspect ratios e as indicated in the plot. The convergence rate becomes independent of e 
for small aspect ratios. 

As an aside, we observe that the remaining error in the approximation 
after the solver has stopped converging, increases with decreasing e. We as- 
cribe this to numerical round off error (all runs are done in double precision). 
The scaling of the condition number with e is consistent with the observed 
loss in accuracy. The smallest used e of 0.0001 still gives an accuracy that is 
acceptable. This remaining error will show up in all further figures. 

In figure [7} the error is plotted against iteration for the vertically lumped 
preconditioner using our additive smoother. We note that the error again 
decays exponentially with iteration number at a rate which is independent 
of e for small aspect ratios. However, one sweep of the additive smoother is 
more expensive than one SOR sweep, and hence it is necessary to compare 
the performance of the two smoothing strategies in terms of computational 
cost as well as number of iterations. 
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Figure 7: Plot showing error (using the oo-norm) against number of iterations, for the 
Conjugate Gradient method applied to the Poisson equation discretised on mesh A, using 
the vertically lumped preconditioner combined with additive smoothing. The mesh has 
been rcscalcd to various different aspect ratios e as indicated in the plot. The convergence 
rate becomes independent of e for small aspect ratios, but the smoother does not improve 
the convergence much for this mesh. 
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In figure [8} the error is plotted against number of iterations for the 
smoothed aggregation preconditioner and the vertically lumped precondi- 
tioner with and without the additive smoother for the matrix obtained from 
Mesh A with e = 0.001. We observe that the vertically lumped precondi- 
tioner converges much faster than the smoothed aggregation preconditioner, 
and that the additive smoother reduces the number of iterations required for 
convergence. The vertically lumped preconditioner has made it feasible to 
solve the pressure Poisson equation on this type of mesh. However, as we 
observe in figure |9] which shows the error plotted against flops, the additive 
smoother is more expensive than SOR since it involves several applications 
of SOR at different levels. We observe that for Mesh A, the time to con- 
verge is approximately the same with the SOR smoother or with the additive 
smoother. The vertically lumped preconditioner produces an approximation 
to the vertically-independent (barotropic) component of the solution with 
very small eigenvalues and it is the job of the smoothers to approximate the 
vertically- varying (baroclinic) component. These results show that for Mesh 
A, which has roughly isotropic tetrahedra before rescaling to small aspect ra- 
tio, the SOR smoother is reasonably effective in approximate the baroclinic 
components. 

Next we present results for Mesh B which is a multiscale mesh, as illus- 



trated in figure 10 In figure 11, we plot the convergence of the smoothed 
aggregation multigrid method applied to the matrix obtained from mesh B, 
which again shows a convergence plateau which becomes longer as e ^ 0. 
Figures 12 and 13 show the convergence rate in iterations for the vertically 
lumped preconditioner with and without the additive smoother respectively. 
In both cases the convergence rate becomes independent of e for small as- 
pect ratios. In this case the additive smoother is producing a faster decay of 
error as the number of iterations increases, compared to the SOR smoother. 
This suggests that the additive smoother is more effective at approximating 
the baroclinic components of the solution, which have a complex multiscale 
structure. The additive smoother uses an algebraic multigrid cycle applied to 
the baroclinic components, which operates at several scales simultaneously. 

Finally in figure 14, the error is plotted against number of iterations for 
the smoothed aggregation preconditioner and the vertically lumped precon- 
ditioner with and without the additive smoother for the matrix obtained 
from Mesh B with e = 0.001. We observe again that the vertically lumped 
preconditioner converges much faster than the smoothed aggregation pre- 
conditioner. Here the vertically lumped preconditioner does not exhibit a 
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Figure 8: Plot showing error (using the oo-norm) against number of iterations, for the 
Conjugate Gradient method appHed to the Poisson equation discretised on mesh A (e = 
0.001), with various different preconditioners. The continuous line indicates the vertically 
lumped preconditioner with the additive smoother, the dashed line indicates the vertically 
lumped preconditioner without the additive smoother, and the dash-dotted line indicates 
the smoothed aggregation multigrid preconditioner. 
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Figure 9: Plot showing error (using the oo-norm) against number of floating point op- 
erations (flops) counted by PETSc, for the Conjugate Gradient method apphed to the 
Poisson equation discretised on Mesh A (e = 0.001), with various different precondition- 
ers. The continuous line indicates the vertically lumped preconditioner with the additive 
smoother, the dashed line indicates the vertically lumped preconditioner without the ad- 
ditive smoother, and the dash-dotted line indicates the smoothed aggregation multigrid 
preconditioner. 
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Figure 10: Plot showing "cutaway" surface through mesh B, used for benchmarking the 
preconditioner. The mesh has very fine mesh elements at the middle of the domain, so 
that there are a large range of lengthscales in the mesh. 

convergence plateau but does have a slow rate of convergence which we at- 
tribute to the presence of the small eigenvalues associated with small scales 
in the mesh which are not altered by the vertically lumped preconditioner. 
The inclusion of the additive smoother means that the number of iterations 
is dramatically reduced, since the additive smoother is a multigrid precon- 
ditioner which treats all of the scales in the mesh. Despite the added cost 
of the additive smoother, we observe that it results in a much more efficient 
solver. We conclude that the additive smoother should be used when small 
scales are present in the mesh which lead to eigenvalues which are of the 
same size as those associated with eigenvectors corresponding to horizontal 
modes. 

7. Summciry and outlook 

In this paper we discussed the ill conditioning of the linear system ob- 
tained from the finite element approximation of the pressure Poisson equation 
on general vertically unstructured meshes in small aspect ratio domains (such 
as the global oceans). We showed that the condition number scales like 
as e — > (where e is the aspect ratio H/L) in the case in which Neumann 
boundary conditions are set on all surfaces. We also showed that the condi- 
tion number is independent of e when Dirichlet conditions are applied at the 
top surface. This motivated a preconditioner consisting of two stages: in the 
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Figure 11: Plot showing error (using the oo-norm) against number of iterations, for the 
Conjugate Gradient method applied to the Poisson equation discretised on mesh B (shown 



in figure 10 e = 0.001), using the smoothed aggregation multigrid preconditioner. The 
mesh has been rescaled to various different aspect ratios e as indicated in the plot. The 
number of iterations required to converge increases for decreasing e, with a long "plateau" 
for small aspect ratios. 
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Figure 12: Plot showing error (using the oo-norm) against number of iterations, for the 
Conjugate Gradient method applied to the Poisson equation discretised on mesh B (shown 



in figure 10 e = 0.001), using the vertically lumped preconditioner. The mesh has been 
rescaled to various difl^erent aspect ratios e as indicated in the plot. The convergence rate 
becomes independent of e for small aspect ratios. 
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Figure 13: Plot showing error (using the cx^-norm) against number of iterations, for 
the Conjugate Gradient method apphed to the Poisson equation discretised on mesh B 
(shown in figure 10 1 , using the vertically lumped preconditioner combined with additive 
smoothing. The mesh has been rescaled to various different aspect ratios e as indicated 
in the plot. The convergence rate becomes independent of e for small aspect ratios, with 
a substantial improvement over the vertically lumped preconditioner with standard SOR 
smoothing, shown in figure |12[ 
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Figure 14: Plot showing error (using the oo-norm) against number of iterations, for 
the Conjugate Gradient method apphed to the Poisson equation discretised on mesh B 
(shown in figure 10 1, with various different preconditioners. The continuous Hne indicates 
the vertically lumped preconditioner with the additive smoother, the dashed line indicates 
the vertically lumped preconditioner without the additive smoother, and the dash dotted 
line indicates the smoothed aggregation multigrid preconditioner. 
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Figure 15: Plot showing error (using the oo-norm) against number of floating point op- 
erations (flops) counted by PETSc, for the Conjugate Gradient method appHed to the 
Poisson equation discretised on Mesh B (shown in figure 10 1, with various different pre- 
conditioners. The continuous line indicates the vertically lumped preconditioner with the 
additive smoother, the dashed line indicates the vertically lumped preconditioner with- 
out the additive smoother, and the dash dotted line indicates the smoothed aggregation 
multigrid preconditioner. 
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first stage an approximate reduced system for the surface degrees of freedom 
is solved, and in the second stage the solution is reconstructed throughout the 
domain with the approximate surface solution used as a Dirichlet boundary 
condition. The first stage results in a much smaller linear system, and the 
second stage involves a submatrix which has a condition number which is in- 
dependent of e. The reduced system is obtained using an algebraic multigrid 
prolongation operator which approximates the vertical extrapolation oper- 
ator, and the second stage submatrix can be solved using further algebraic 
multigrid stages. Using numerical experiments, we showed that this precon- 
ditioner, when combined with the Conjugate Gradient method, results in a 
solver which has a convergence rate that is independent of the aspect ratio. 
Further, we showed that the additional computational cost of using the ad- 
ditive smoother means that it is only beneficial in truly multiscale meshes. 
Those are meshes that do not just have two entirely different lenght scales, 
the horizontal and the vertical, but a whole range of scales inbetween . This 
strategy will become crucial when solving process study problems consisting 
of small scale dynamics (such as open ocean deep convection, or density over- 
flows) that are embedded in a large scale circulation. We also anticipate that 
the smoother will become important when parallel domain decomposition 
methods are used, where (block) SOR smoothing methods are known to be 
less effective. The methods described in this paper have been implemented 
in ICOM where they will be used to investigate adaptive unstructured large 
scale ocean modelling. 
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